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Weakly Interacting Massive Particles (WIMPs) are one of the main candidates for the dark matter 
in the Universe. If these particles make up the dark matter, then they can be captured by the Sun or 
the Earth, sink to the respective cores, annihilate, and produce neutrinos. Thus, these neutrinos can 
be a striking dark matter signature at neutrino telescopes looking towards the Sun and/or the Earth. 
Here, we improve previous analyses on computing the neutrino yields from WIMP annihilations in 
several respects. We include neutrino oscillations in a full three-flavor framework as well as all effects 
from neutrino interactions on the way through the Sun (absorption, energy loss, and regeneration 
from tau decays). In addition, we study the effects of non-zero values of the mixing angle #13 as well 
as the normal and inverted neutrino mass hierarchies. Our study is performed in an event-based 
setting which makes these results very useful both for theoretical analyses and for building a neutrino 
telescope Monte Carlo code. All our results for the neutrino yields, as well as our Monte Carlo code, 
are publicly available. We find that the yield of muon-type neutrinos from WIMP annihilations in 
the Sun is enhanced or suppressed, depending on the dominant WIMP annihilation channel. This 
effect is due to an effective flavor mixing caused by neutrino oscillations. For WIMP annihilations 
inside the Earth, the distance from source to detector is too small to allow for any significant amount 
of oscillations at the neutrino energies relevant for neutrino telescopes. 

PACS numbers: 14.60.Pq, 95.85. Ry, 95.35.+d 



I. INTRODUCTION 

Weakly Interacting Massive Particles (WIMPs) are interesting (and perhaps the most plausible) candidates for the 
dark matter in the Universe. In particular, neutralinos (that arise naturally in supersymmetric extensions of the 
Standard Model) are promising WIMP candidates and can be captured by e.g. the Sun or the Earth as the solar 
system moves in the Milky Way halo [3, [3, Once captured, they will sink to the core and annihilate into several 
different annihilation channels (leptons, quarks, gauge bosons, or Higgs bosons). When the annihilation products 
decay, neutrinos of different flavors are produced (in the case of e.g. Kaluza-Klein dark matter, neutrinos can also 
be produced directly in the annihilation processes [1]). During the propagation out of the Sun or the Earth, these 
neutrinos can interact and undergo neutrino oscillations. Here, we will improve on earlier estimates of the neutrino 
yields from WIMP annihilations in the Sun and the Earth by performing a very detailed analysis including both 
neutrino oscillations and interactions. 

Neutrino oscillations will be included in a full three-flavor framework including both vacuum and matter oscillations 
all the way from the production region to a possible detector. Interactions are only relevant for neutrinos propagating 
out of the Sun and we will include absorption via charged-current (CC) interactions, energy losses via neutral-current 
(NC) interactions, and regeneration of neutrinos from tau decay (after CC interactions). Both neutrino oscillations 
and interactions will be treated simultaneously in a consistent way. 

We will perform this analysis in an event-based framework, meaning that we will follow each neutrino from pro- 
duction to detection. The advantage of this approach is that, at the same time as we obtain theoretical estimates of 
the neutrino yields at a neutrino telescope, we will also produce events that can be used in Monte Carlo simulations 
by the neutrino telescope community. In addition, our code, which is written in the programming language Fortran, 
is publicly available and can be implemented along with a neutrino telescope Monte Carlo, which means that such a 
simulation is not restricted to our finite list of parameter choices. 
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A recent analysis by Cirelli et al. [5| has been performed along the same lines as our work. The main difference 
between their work and our work is that they focus on the neutrino flavor distributions based on a density matrix 
formalism and do not follow each neutrino individually from the production point to the detector. We will make 
comparisons with their results later in our paper. In the final stages of preparation of this paper, we became aware 
of the studies by Lehnert and Weiler 0] as well as Barger et al. 0]. In the former, the neutrino flavor flux ratios from 
WIMP annihilations in the Sun in the absence of neutrino interactions are presented (which is valid up to neutrino 
energies of the order of 10 GeV, above which neutrino interactions affect the fluxes by more than a few percent). We 
will find that our low-energy results agree very well with those of Lehnert and Weiler, while our high-energy results 
differ due to our inclusion of neutrino interaction effects. In the latter, the approach is similar to that of Cirelli et al, 
but also includes the dependence on the spin of the WIMPs and mainly discusses annihilation channels where equal 
amounts of all neutrino flavors are produced, and thus, the neutrino oscillation effects are suppressed. Furthermore, 
our study also includes neutrinos created by WIMP annihilations in the Earth. In addition, Crotty @ has earlier 
published the results of using a similar Monte Carlo code when treating the annihilation of supermassive (10 8 GeV 
< m x < 10 16 GeV), strongly interacting dark matter particles. 

This paper is organized as follows. In Sec. [TTJ we will review how WIMPs are captured in the Sun or the Earth and 
how annihilations are treated. Next, in Sec. IIII1 we will describe our framework for including neutrino interactions 
and oscillations into our Monte Carlo code. Then, in sections IIVIIVI1 we will discuss how we perform the neutrino 
propagation out of the Sun, from the Sun to the Earth, and through the Earth to an actual detector. The intermediate 
results for the neutrino yields of different flavors at the surface of the Sun, at a distance of 1 AU, and at the actual 
detector will be presented. In Sec. lVIIl we will investigate the impact of non-zero #13 as well as that of having normal 
or inverted neutrino mass hierarchy, and in Sec. IVIlTl we will consider WIMP annihilations in the Earth. Finally, in 
Sec. IIX1 we will summarize our results and present our conclusions. In addition, in the appendix, we will show the 
equivalence of our Monte Carlo simulation and the density matrix formalism (used in e.g. Ref. [5j). 



II. ANNIHILATION OF WIMPS IN THE SUN/EARTH 

WIMPs can be captured in the Sun |l[ or the Earth [2j when they move through the dark matter halo in the 
Milky Way. Once the WIMPs have been captured, they will sink to the core of the Sun or the Earth, where they 
can annihilate and produce leptons, quarks, gauge and Higgs bosons. These will hadronize, decay, and eventually, 
produce neutrinos. We therefore follow the calculations performed in Refs. [H, [To| . but update them to be more 
general and include more annihilation channels. We use Pythia 6.400 [11] to simulate the hadronization and decay 
of the annihilation products and collect the neutrinos and antineutrinos produced. The annihilation channels that 
we simulate are the following: cc, bb, tt, t + t~, W + W~ , Z°Z°, gg, v e v e , v^v^, and v T v T [HI], where the last three 
channels are only relevant for WIMP candidates that can annihilate directly to neutrinos (e.g. Kaluza-Klein dark 
matter) . The lighter charged leptons are not of importance, since electrons are stable and muons are stopped before 
they have a chance to decay and produce (high-energy) neutrinos. For the bb channel, we need to be especially careful, 
since these hadronize and produce B mesons, which, for the case of the Sun, interact before they decay. We include 
these interactions in an approximate fashion by performing the Pythia simulations as if in free space, and later, we 
rescale the energy of the resulting neutrinos by estimating the energy loss of the B mesons due to their interactions. 
This follows the calculation in Ref. [12j, but we use newer estimates of the B meson interaction cross-sections (as 
given in Ref. [IH). For the solar density, we use the Standard Solar Model [3]. Annihilation channels including Higgs 
bosons are not simulated directly, since we can easily calculate the neutrino yield from the decay of the Higgs bosons 
given the Higgs masses and decay branching rates. For a specific WIMP model, we later perform this step separately, 
once the Higgs properties are known. In this calculation, we let the Higgs bosons decay in flight, properly integrating 
over the decay angles and performing the appropriate Lorentz boost to obtain the neutrino yields in the rest frame 
of the Sun or the Earth. We perform simulations for the following WIMP masses: 10, 25, 50, 80.3, 91.2, 100, 150, 
176, 200, 250, 350, 500, 750, 1000, 1500, 2000, 3000, 5000, and 10000 GeV. For each mass and annihilation channel, 
we perform 2.5 x 10 6 annihilations and collect all neutrinos and antineutrinos, keeping track of the flavors separately. 
In the next section, we will describe our formalism for neutrino interactions and oscillations, and in the succeeding 
sections, we will discuss how we propagate the neutrinos from the core of the Sun and/or Earth to the detector. 

In Fig. [TJ we present our results for the neutrino yields at creation in some of the different annihilation channels. 
Comparing this to Fig. 2 in the published version of Ref. [5j, we find that our results for the W + W~ annihilation 
channel is about a factor of two lower and that there is a discrepancy in the tt channel, where we also have a lower 
yield. We also see a double-peaked yield for the tt channel compared to the broader single-peaked yield in Ref. [f|. 
In addition, our high-energy peak goes up to higher energies than their results. This two-peak structure is expected 
as the high-energy peak comes from the prompt decay of W bosons and the low-energy peak comes from quark jets. 

We have been in contact with the authors of Ref. [5| regarding these differences, and as a consequence of this, they 
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FIG. 1: The neutrino yields for electron and tau neutrinos as functions of z = E v /m x for six different WIMP annihilation 
channels at production in the center of the Sun and the Earth. Note that the muon neutrino yields are the same as the electron 
neutrino yields and are therefore not shown separately. 



have found two errors in their code. The first one was a factor of two too high yields for the W + W~ channel (which 
affected the yield from W + W~ channel directly and the yield from tt indirectly). The other one was an error in the 
Lorentz boost of the W boson resulting from top decay. After the authors of Ref. || corrected these errors, our yields 
now agree reasonably well [44|. There is still a small difference in the tt channel, but this difference is mostly due to 
our inclusion of final state radiation of gluons. 



III. NEUTRINO INTERACTIONS AND OSCILLATIONS 



On their way from creation in the center of the Sun or the Earth to measurement in the detector, neutrinos can 
undergo both neutrino-nucleon interactions and neutrino oscillations. These are treated according to the following 
two prescriptions. 



A. Neutrino interactions 



Neutrinos can undergo both CC and NC interactions with nucleons during their passage through a medium. In 
a CC interaction, we obtain a charged lepton, whereas in a NC interaction, we obtain a neutrino with degraded 
energy. We calculate the total and differential cross-sections for these interactions, as both are needed to determine 
the absolute interaction rates and the energy loss of the lepton. We need the latter to take into account energy losses 
of neutrinos in the Sun and regeneration of neutrinos from tau decay, as well as for calculating the scattering process 
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near or in the neutrino detector. We use the convenient expressions for the (anti)neutrino-nucleon cross-sections that 
can be found in Ref. [l5|. These expressions include the effect of the tau mass on the (anti)neutrino-nucleon CC 
cross-sections. 

In order to calculate the cross-sections, we need to specify the parton distribution functions (PDFs). We use the 
CTEQ6 [l6| PDFs, or more specifically, the CTEQ6-DIS PDFs. For the other physical parameters, we choose the 
following values 



GeV~ 



M w 


= 80.41 GeV, 


M z 


= 91.188 GeV, 


G F 


= 1.16639 x 1(T 5 


M p 


= 0.938272 GeV, 


M n 


= 0.939566 GeV, 


■ &w 


= 0.23124. 



The resulting cross-sections agree well with other calculations in this energy range, see e.g. Ref. [17j (however, only 
cross-sections on isoscalar targets are given in that work) . 

The neutrino-nucleon CC cross-sections have been measured by several experiments (see Ref. [l8[ for a summary). 
In the energy range up to E v = 350 GeV, the average cross-sections are [18j 

— = (0.677 ± 0.014) x 10~ 38 cm 2 /GeV, 
E v 

^ = (0.334 ± 0.008) x 10" 38 cm 2 /GeV. 
As a comparison, our calculated cross-sections at 100 GeV are 



= 0.684 x 10 - 38 cm 2 /GeV, 



E v = 100 GeV 
1 _f2i_ = 0.330 X 10- cmVGeV, 

i.e., very close to the measured cross-sections. 

Note that, in order to gain computational speed, we do not use the full expressions for the total cross-sections 
directly, instead we calculate the different cross-sections for a range of energies from 1 GeV to 10 12 GeV and interpolate 
these results. The interpolation errors introduced are below 1 %, but the speed increase is dramatic. 

In order to take care of each neutrino-nucleon interaction, we have developed a Monte Carlo code that simulates the 
interaction and returns the energy and angles of the final state lepton as well as the hadronic shower. The neutrino- 
interaction and Monte Carlo codes are available for download at the website in Ref. 19], where more technical 
information is also available. 

Our differential spectra differ somewhat from those obtained with e.g. Pythia. The main reason is that Pythia 
does not allow us to use the full phase-space (i.e., going down to sufficiently low momentum transfer Q 2 ). This is not 
surprising, since this is not a process for which Pythia is optimized. We can also observe this effect looking at the 
total cross-sections, which Pythia underestimates at low energies (less than 100 GeV), again due to the fact that the 
phase-space is not complete. 



B. Neutrino oscillations 



We want to describe the yields of neutrinos at the Earth as accurately as possible. In order to accomplish this 
task, we need to include the effects of neutrino oscillations, which have been observed in experiments (20l . l2ll I22I I23I 
HH, |H, |2f|, H3, [H, HH . For the high-energy neutrinos that we are considering, the adiabatic propagation of neutrino 
matter eigenstates inside the Sun, which is used for ordinary solar neutrinos, is no longer valid. The reason is mainly 
that the neutrinos are produced at densities which are larger than the high-energy resonance density, corresponding 
to the resonance of the small vacuum mixing angle #13. In addition, for certain setups (most notably for high-energy 
neutrinos), it may not be a good idea to use the approximation of decoherent mass eigenstates arriving at the Earth, 
which is why we store these neutrinos on an event basis rather than as energy spectra for the different mass eigenstates. 

The theory of neutrino oscillations is based on the assumption that the neutrino flavor eigenstates are not equivalent 
to the neutrino mass eigenstates, but rather that they are linear combinations of each other. This can be written as 
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where U is the leptonic mixing matrix, \v a ) is a neutrino flavor eigenstate, and \v a ) (a = 1, 2, 3) are the neutrino mass 
eigenstates with definite masses m a . The evolution of the neutrino state v(t) — (u e (t), ^u(t), v T (t)) T is given by 

v{t) = S(t)u(Q), (2) 

where the evolution operator S(t) depends on the distance and medium traversed. For three-flavor neutrinos prop- 
agating in matter of constant density, the evolution operator is given by S(t) = exp(— iHt), where H is the total 
Hamiltonian including both the vacuum and the Mikheyev-Smirnov-Wolfenstein (MSW) potential terms, i.e., 

H = -3-J7diag(0, Am^, Amfi)^ + di&g(V2G F N e ,0,0). (3) 

Here E is the neutrino energy, Am 2 ^ = mf — m| are the neutrino mass squared differences, Gf is the Fermi coupling 
constant, and N e is the electron number density. By the use of the Cay ley-Hamilton formalism, it is possible to write 
the evolution matrix as (30j 

{)= *h — ' (4) 

where <fi is an overall phase factor which does not affect neutrino oscillations, / is the 3x3 identity matrix, T = 
H — tr(P)//3 is the traceless part of the Hamiltonian, c\ = — tr(T 2 )/2, and A a are the eigenvalues of T [which can 
be determined from ci and det(T)]. The probability of an initial neutrino of flavor v a to be in the flavor eigenstate 
vp at time t is then given by P a p{t) = \S0 a (t)\ 2 . 

When propagating neutrinos through matter of varying density, we can approximate the electron number density 
profile by a large number of layers with constant electron number density (the appropriate number of layers is 
determined by the rate of change in the electron number density as well as the neutrino oscillation lengths). If we let 
a neutrino pass through k layers and label the evolution operator of layer i by Si, then the total evolution operator 
S is given by [3l| 

S = SfeSfc-i . . . S 2 Si. (5) 

Neutrino oscillations in vacuum can also be treated with this method, we simply use one layer with the appropriate 
length and set the electron number density to zero. We will later calculate S with these methods using accurate 
matter density profiles for both the Sun and the Earth. Note that the calculation of S in Eq. ([5]) is particularly simple 
given the Cayley-Hamilton formalism presented in Eq. 

The neutrino oscillation probabilities only depend on the leptonic mixing matrix and the mass squared differences. 
The standard parameterization of the leptonic mixing matrix is given by [32| 

C13C12 C13S12 

U = \ -S12C23 - Ci 2 S23Sl3e 1<5 C12C23 - Si2S23Sl3e I<5 S23C13 I , (6) 
S12S23 - Cl2C23Sl3e" 5 -C12S23 - Sl2C23Sl3e 115 

where cy = cos(#y ), Sij = sin(#y ), and S is the CP- violating phase. The present phenomenological constraints on 
the neutrino oscillation parameters are given by global fits to the different neutrino oscillation experiments as [33[ 

012 = 33.2° ±4.9°, 

013 < 12.5°, 

023 = 45.0° ± 10.6°, 
Se [0,2tt), 

Am^ = (8.]+S;§) x lCT 5 eV 2 , 
|Am 2 1 | = (2.2l^)xlO- 3 eV 2 , 

where the presented values are the best-fits and the 3er (99.7 % confidence level) ranges. The leptonic mixing angle 0i 2 
and the small mass squared difference Ato^i are mainly constrained by solar and long-baseline reactor experiments, 
while the leptonic mixing angle 023 and the large mass squared difference Ato 2 ^ are mainly constrained by experiments 
studying atmospheric and accelerator neutrinos. The upper bound on the leptonic mixing angle 0i3 is currently given 
by short-baseline reactor experiments, while the CP-violating phase 5 is completely unknown. 

In our simulations, we have focused on the best-fit values for 0i 2 , 023, Am 2 i, |Am|i|, and put the CP- violating 
phase S equal to zero. The different combinations of the following scenarios were implemented: 
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• Three different values of #13 (#13 = 0, 5°, 10° corresponding to sin 2 26\^ = 0, 0.03, 0.12, respectively) 

• Normal or inverted mass hierarchy (i.e., Am^ positive or negative) 

In addition, we have implemented the case when there are no neutrino oscillations by setting fty = and Amf 3 - = 0. 
The reason is to study the effect of neutrino oscillations separately as well as being able to compare our non-oscillation 
results with previous works that have not taken neutrino oscillations into account. For the remainder of this paper, 
the examples shown are for $13 = and normal neutrino mass hierarchy unless stated otherwise (c/., Sec. IVII|) . 
Thus, our standard set of neutrino oscillation parameters (std. osc.) is #12 = 33.2°, #13 = 0, #23 = 45°, 6 = 0, 
Aml x = 8.1 x 10~ 5 eV 2 , and Am§ x = 2.2 x 10" 3 eV 2 . 

IV. PROPAGATION OUT OF THE SUN 

We are now ready to perform the steps necessary to propagate the neutrinos produced at the center of the Sun or 
the Earth calculated in Sec. [TTJ We start with the Sun, as it is more involved, and take care of the Earth later. 

For the Sun, we divide the propagation to the detector into several steps, starting with the propagation from the 
center of the Sun to its surface. In this process, we take into account both neutrino interactions (CC and NC) as well 
as neutrino oscillations. 



A. Solar model 

Let us first describe the solar model that we have used. We need a good knowledge of the Sun's interior in order to 
know both the matter effects on neutrino oscillations and the interaction probabilities. We have used the Standard 
Solar Model [l4[ , which gives the density of the dominant elements and electrons as a function of the radial distance 
from the center of the Sun. We use the element densities to calculate the density of protons and neutrons as a function 
of radius, which is needed for the neutrino interaction probabilities, and use the electron densities to calculate the 
matter potential needed for the neutrino oscillation matter effects. 

B. Neutrino interactions and oscillations 

Since neutrino interactions and oscillations can occur simultaneously, we need to take both into account at the 
same time. However, the neutrino-nucleon cross-sections are practically flavor independent, which means that this is 
particularly easy to perform. (This is true as long as the tau mass can be neglected. We include the effects of the 
tau mass later on when we treat the tau neutrino CC interactions.) For each neutrino, we first calculate a creation 
point. Since the WIMPs are in thermal contact with the Sun's core, they will be distributed according to a Gaussian 
distribution 

n(r) = n(0)e~ r2 /( 2r x) (7) 

with 

( 3kT \ 1/2 

where k is the Boltzmann constant, T is the core temperature in the Sun, G is the gravitational constant, p is the 
core density, and m x is the WIMP mass. For the Sun, annihilations are concentrated to within about 1 % of the 
solar radius, which means that the effect is small. We choose a creation point and calculate the path length from 
this creation point to the surface of the Sun. Due to the concentration of the annihilations to the very core, we can 
safely approximate this path as being purely radial to simplify our calculations. After this stage, we know how much 
material the neutrino has to traverse in order to reach the solar surface. We then compute the neutrino-nucleon 
cross-section and use it to randomize the point of interaction for each neutrino (or if the neutrino escapes the Sun 
without interacting). Depending on the type of interaction (CC or NC), we simulate the interaction results. 

In the case of a CC interaction, we use the neutrino evolution operator method described in Sec. MI Bl to determine 
the probabilities of the neutrino being in the different flavor eigenstates at the time of interaction. The actual 
interacting neutrino flavor is then randomized according to these probabilities. If the neutrino at this interaction 
point is a tau (anti)neutrino, we simulate the neutrino interaction (properly including the tau mass suppression of the 
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FIG. 2: The neutrino yields as a function of z = E v /m x at creation in the center of the Sun (left panel) and at the surface of 
the Sun (right panel) for annihilation of 250 GeV WIMPs into r + r _ . The secondaries in the right-hand panel are neutrinos 
that come from r decay after a v T jv T CC interaction. 



tau (anti)neutrino-nucleon cross-section [45l|). We then obtain a charged tau lepton, whose decay is simulated with 
Pythia. We collect the neutrinos produced in this decay and proceed with their propagation out of the Sun. If the 
neutrino is an electron or muon neutrino, then the charged lepton produced in the CC interaction is either stable or 
will be stopped before it has had time to decay and the neutrino is considered to be absorbed. 

The main effect of a NC interaction is to change the energy of the neutrino. Since the NC vertex is flavor-blind, 
it does not destroy the relative phases among different neutrino flavor eigenstates and does not force the neutrino 
to interact as a particular flavor (see e.g. Ref. (34[ for a discussion on the oscillations of neutrinos from Z°-decay, 
the argument can be readily extended to NC interactions of neutrinos). After an interaction, we then repeat this 
procedure until the neutrino has reached the surface of the Sun. In App. [A] it is shown that the Monte Carlo method 
described above is statistically equivalent to the density matrix formalism used in Ref. 

For the neutrino oscillations between creation and interaction points, we use the method described in Sec. IIIIBI 
This means that we divide the Sun, with its varying matter density, into layers of constant density. We have chosen 
to let each layer have a width of 0.3 % of the radius of the Sun. This gives a total error on the neutrino state after 
propagation of much less than 1 %, which is a reasonable compromise between accuracy and speed. 



C. Results for neutrino yields at the surface of the Sun 

We are now ready to present some results at the solar surface in order to show the effects of propagation through 
the Sun. Since most annihilation channels produce similar amounts of neutrinos of different flavors, we focus on 
annihilation to r + r _ , which produces far more tau neutrinos than electron and muon neutrinos, in order to emphasize 
the effects of neutrino oscillations. In Fig. [3J we show the yields of neutrinos from annihilation of 250 GeV WIMPs 
to t + t~ both at the center and at the surface of the Sun. This figure is derived for the standard set of neutrino 
oscillation parameters as described earlier. We can clearly observe two main effects in this figure. 

The first effect is that oscillations effectively mix muon and tau neutrinos on their way out of the Sun, whereas 
the electron neutrinos remain essentially unmixed. This is due to the electron neutrino being equivalent to one of 
the matter eigenstates when propagating in a medium of high electron number density (which is true for the most 
part of the neutrino propagation inside the Sun). The almost complete mixing of muon and tau neutrinos is due to 
the effective mixing angle in the remaining two-flavor system being very close to 623 [IH, which is maximal in our 
simulations, while the dependence on 612 is suppressed by the ratio a = Amli/Am^. 

The second effect is the pile-up of events at low energies. Some of these events come from secondary decays of 



8 



10 



M. Blennow, J. Edsjd and T. Ohlsson, 2007 




10 



-2 



10 



-3 



10 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



z = E v / m x 



FIG. 3: The muon neutrino and antineutrino yields as a function of z = E v /m x at creation and at the surface of the Sun for 
annihilation of 250 GeV WIMPs into W + W~ . In this figure, neutrino oscillations have been turned off in order to show the 
effect of neutrino interactions. 

that have been created in v T jv T CC interactions and some of these come from energy losses due to NC interactions. 
The interactions also create a loss at higher energies (both due to absorption via CC interactions of electron and muon 
neutrinos and due to energy losses). This is not fully visible in this figure, but in Fig. [31 we can observe this effect 
more clearly for the W + W~ channel with neutrino oscillations turned off. The loss of neutrinos at higher energies 
due to these interactions is clearly visible. 



When the neutrinos have escaped from the Sun, we use the neutrino evolution operator method described in 
Sec. MI Bl to propagate the neutrinos from the surface of the Sun to a distance of I AU from the Sun. In the treatment 
of the oscillations of ordinary solar neutrinos (produced by the thermonuclear reactions sustaining the Sun), it is 
common practice to assume that the neutrinos arrive at the Earth as mass eigenstates. The main arguments for 
making this assumption are the following: 

• The coherence length due to separation of wave packets for the neutrino mass eigenstates is much shorter than 



• The change in the neutrino baseline due to the diameter of the Earth and the eccentricity of the Earth's orbit 
is comparable to or larger than the oscillation length and also gives an effective average of the oscillations. 

• The finite energy resolution of detectors results in an effective averaging of the resulting fast oscillations. 

The first of these assumptions is no longer valid when we study the high-energy neutrinos produced by WIMP 
annihilations. Thus, we save the neutrinos arriving at 1 AU with both energy and relative phase information among 
the neutrino flavor eigenstates (in order to make it possible to propagate the neutrinos further), rather than making 
energy spectra for the different flavors. Depending on the actual neutrino flux, neutrino telescopes may be sensitive to 
its temporal variation, which essentially would render the second assumption invalid (we will come back to this issue 
later). However, the energy resolution of a typical neutrino telescope is still not very high, making it plausible that the 
third assumption could still hold (in which case we could just compute the energy spectra for the different neutrino 
mass eigenstates arriving at the Earth). However, in some dark matter scenarios, such as Kaluza-Klein dark matter, 
the WIMPs may annihilate directly into pairs of neutrinos and antineutrinos, resulting in a sharp monochromatic 
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FIG. 4: The neutrino yields as a function of z = E u /m x at 1 AU for annihilation in the Sun of 250 GeV WIMPs into t + t~ 
both with (std. osc.) and without neutrino oscillations. The main effects of oscillations are easily seen, i.e., muon and tau 
neutrinos are fully mixed, and electron neutrinos are partially mixed with the other flavors. 

peak in the spectrum. The smudging of the energy spectrum from the poor energy resolution of the detector does not 
change the neutrino oscillation probabilities of such monochromatic neutrinos, and thus, we keep the most general 
setup in order for our method to be independent of the type of dark matter that is studied. The results after this 
stage are stored on an event-by-event basis complete with neutrino energies, amplitudes and phases. 

A. Results for neutrino yields at 1 AU 

The vacuum neutrino oscillations from the Sun to the detector force the electron neutrinos to mix with the muon 
and tau neutrinos. In Fig. 2J we show the yields at 1 AU for 250 GeV WIMPs annihilating into t + t~ both with 
and without neutrino oscillations. Compared to the results at the surface of the Sun (c/., Fig. [2]), we note that also 
the electron neutrinos are effectively mixed with the muon and tau neutrinos. This effect is solely due to vacuum 
neutrino oscillations governed by the small mass squared difference Amji (since 6*13 = 0). As can be seen, the mixing 
is large, but not complete, which is due to the leptonic mixing angle #12 being large but not maximal. Furthermore, 
the yields of muon and tau neutrinos now essentially coincide due to the fast oscillations governed by the neutrino 
oscillation parameters Am^ and 6*23 (because of the L/E dependence of the neutrino oscillation phase, the phase 
is more sensitive to the neutrino energy E at a distance of 1 AU than at the surface of the Sun). When comparing 
the results with and without oscillations, we can observe that the main impact of the oscillations is to weaken, or 
even erase, asymmetries in the initial yields through lepton mixing. As a result of this, any WIMP candidate, which 
produces an excessive amount of tau neutrinos, will have its resulting muon neutrino yield at a distance of 1 AU 
increased (in this case by a factor of 3-4). Thus, it will be easier to perform indirect detection of such a WIMP 
candidate than what may be naively expected if not including neutrino oscillations. On the other hand, any WIMP 
candidate, which produces an excessive amount of muon neutrinos, will have its muon neutrino yield at a distance 
of 1 AU decreased and will therefore be harder to detect. Furthermore, it is interesting to note that, while electron 
and muon neutrinos are produced in equal amounts, neutrino oscillations affect the yields in such a way as to have 
an equal amount of muon and tau neutrinos in the later stages of evolution. 
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FIG. 5: The time-averaged neutrino yields as a function of z = E v /m x at the South Pole for annihilation in the Sun of 250 GeV 
WIMPs into t + t~ both with (std. osc.) and without neutrino oscillations. 

VI. PROPAGATION TO AND THROUGH THE EARTH TO AN ACTUAL DETECTOR 

In the preceding section, we propagated the neutrinos to a distance of 1 AU from the Sun and we are now ready to 
proceed with the final propagation to the Earth and through the Earth to an actual detector. The reason to separate 
the propagation to 1 AU and the final propagation to the detector is mainly practical. The propagation to 1 AU is 
neither time nor detector dependent, but the final propagation to the actual position of the Earth and the detector 
is both time and detector location dependent. Hence, we have split the propagation in this way to make it easier to 
use the results of the previous section for Monte Carlo simulations for actual detector responses. 

In order to include these effects, we use a simple model of the Earth's motion. This model contains a few ap- 
proximations, but still incorporates the interesting physics effects. We include the eccentricity of the Earth's orbit 
and the daily rotation around its own axis, but we make the simplifying assumption that both perihelion and the 
winter solstice occurs at New Year (both are off by about one week, but this approximation has no large effect on 
our results). Within this model, we choose a time (fraction of the year since New Year) of each event and for this 
time we calculate the Earth's actual distance from the Sun. We then use this information to propagate the neutrinos 
from 1 AU to the actual position of the Earth (in case the Earth is closer than 1 AU, we just apply the inverse of the 
propagation operator). In future releases of the WimpSim code, we may include a more sophisticated astronomical 
calculation. This has no effect on the results of this paper, since the errors introduced by these approximations are 
not observable. 

Furthermore, we have to determine the path length traversed through the Earth. We perform this by specifying the 
latitude of the detector and calculate the orientation of the Earth at the time of the event. We use this knowledge to 
compute the path the neutrino will take through the Earth. The propagation through the Earth is done in the same 
way as through the Sun, i.e., we apply the neutrino evolution operator method described in Sec. IIII Bl where we use 
the Earth's matter density profile as given in Ref. [3(| (from which we calculate the electron density as a function 
of radius). We then propagate the neutrino in steps of 0.3 % of the radius of the Earth, which turns out to be a 
good compromise between accuracy and speed (the accuracy is far better than 1 %). Even in this case, we make an 
astronomical calculation in order to simplify our simulation, which is that the detector is pointing towards the Sun 
and January 1 at 00:00. Again, this approximation has no effect on the results of our calculation. 

Once we have propagated the neutrinos to the detector, we let them interact and again simulate their CC and/or 
NC interactions to produce hadronic showers and/or charged leptons. Both the neutrino yields at the detector as well 
as the hadronic showers and neutrino-induced lepton yields can be calculated. In case a charged muon is produced, we 
further let it traverse the detector medium where it can undergo both energy losses and multiple Coulomb-scattering 
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to obtain the muon flux. In order to calculate the neutrino interactions and energy losses of muons properly, we need 
to specify what kind of medium the detector and its surroundings are made of. We have included both ice/water and 
rock as possible mediums. We do not show any results for the hadronic showers or charged leptons at the detector, 
since they contain no new information regarding the neutrino oscillation effects. Nevertheless, all these results are 
available at the website in Ref. [1 91 ] - 

A. Results for neutrino yields in an actual detector 

In order to be more specific, let us assume that the detector is at —90° latitude (e.g. the IceCube detector [37j at 
the South Pole) and that the viewing period is from spring equinox to autumn equinox, i.e., when the Sun is below 
the horizon for a detector at the South Pole. For the same channel as before, t + t~ from 250 GeV WIMPs, we show in 
Fig. the neutrino yields at this detector, averaged over this viewing period. The main effect is that the oscillations 
that were present in previous figures are now almost completely washed out. It turns out that the matter effects from 
the passage through the Earth are negligible and the main reason for this "averaging" comes from the eccentricity of 
the Earth's orbit, i.e., the variation of the Sun-Earth distance effectively causes a wash-out of the oscillation pattern 
seen at 1 AU. Of course, this is an artifact of our choice of averaging the yield over the viewing period between the 
spring and autumn equinoxes in Fig. [5] In addition, neutrino telescopes have a rather poor energy resolution (50 % or 
worse in this energy range), and in practice, the oscillations in energy seen at a specific time (and Earth-Sun distance) 
are not observable anyway. 

However, for a source of monochromatic neutrinos (e.g. annihilation directly to vT>), this effect could in principle be 
observable, if the number of events is large enough. This would then be seen as a variation of the event rate over time 
as the Earth would move in and out of the oscillation phases. Nevertheless, for the leading dark matter candidates, 
neutralinos in supersymmetric models and Kaluza-Klein dark matter, this signature would not be observable. The 
reason is that neutralinos do not annihilate to vv directly, whereas for Kaluza-Klein dark matter that do annihilate 
to vv, they annihilate to all flavors with equal branching ratios, and then, neutrino oscillations will not affect the 
yields anyway. 

VII. NON-ZERO 6>i3 AND COMPARISON BETWEEN NORMAL AND INVERTED MASS 

HIERARCHY 

In this section, we investigate the impact of a non-zero leptonic mixing angle #13 as well as the differences between 
normal and inverted mass hierarchies on the neutrino yields from WIMP annihilations. We also discuss the differences 
that appear when considering antineutrinos. We present the neutrino yields at the surface of the Sun as well as at 
the distance of 1 AU. 

In Fig. [6l we plot the neutrino yields for #13 = as well as for #13 = 10° and normal mass hierarchy at the surface 
of the Sun. We observe that, in the case of a non-zero #13, the yield of electron neutrinos is enhanced compared with 
the yields for muon and tau neutrinos, which are accordingly suppressed. In fact, for #13 = 10°, the averages of the 
yields for different neutrino flavors are comparable in magnitude, especially for low neutrino energies (note that we 
plot the combined muon and tau neutrino yields). The conversion of muon and tau neutrinos into electron neutrinos 
for non-zero #13 is due to adiabatic MSW flavor conversions taking a larger role at the high MSW resonance. For 
#13 = 0, adiabaticity (and thus also the flavor conversion) is lost. In the case of inverted mass hierarchy, the yields for 
#13 = and #13 = 10° are qualitatively equal to the yields for normal mass hierarchy and #13 = 0. This means that 
a non-zero value of the mixing angle #13 has no effect for inverted mass hierarchy. In the case of antineutrinos, the 
situation is the opposite, the yields for normal mass hierarchy are qualitatively the same, whereas there is a difference 
between the yields for inverted mass hierarchy. In this case, the reason is that considering antineutrinos instead of 
neutrinos the replacement V — > —V has to be performed, and thus, there is a resonance for antineutrinos only if 
Am^ < (inverted mass hierarchy). Naturally, the antineutrino yields for normal and inverted mass hierarchies are 
qualitatively the same for #13 = 0. 

In Fig. [7J we plot the neutrino yields corresponding to those in Fig. O but instead at the distance of 1 AU. Again, 
we observe that in the case of normal mass hierarchy, the neutrino yields for #13 = and #13 = 10° arc different from 
each other. In the case of #13 = 0, the electron neutrino yield is generally lower than the yields for muon and tau 
neutrinos, whereas for #13 = 10°, the electron neutrino yield is generally higher than the yields for muon and tau 
neutrinos. As in the case at the surface of the Sun and for low neutrino energies, the averages of the yields of different 
neutrino flavors are comparable in magnitude. However, in the case of inverted mass hierarchy, the neutrino yields 
are essentially independent of the mixing angle #13, and the yields are similar to those in the case of normal mass 
hierarchy and #13 = 0. In the case of antineutrinos, the discussion of the yields at the distance 1 AU is similar to the 
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FIG. 6: The neutrino yields as a function of z — E v /m x for normal mass hierarchy at the surface of the Sun for annihilation 
of 250 GeV WIMPs into r+r" for both 0i 3 = (std. osc.) and Ois = 10°. 
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FIG. 7: The neutrino yields as a function of z = E v /m x for normal mass hierarchy at the distance of 1 AU for annihilation of 
250 GeV WIMPs into r+r" for both 0ia = (std. osc.) and 6ia = 10°. 



discussion for the yields at the surface of the Sun, which means that effects on the yields appear for 6*13 = 10° and 
inverted mass hierarchy. In all other cases, the yields are qualitatively equal. 

We conclude that a non-zero leptonic mixing angle #13 and the neutrino mass hierarchy may influence the neutrino 
yields from WIMP annihilations. Therefore, neutrino yields from WIMP annihilations could, in principle, be used 
as a tool in determining both a non-zero value of the mixing angle #13 as well as the mass hierarchy of the neutrino 
mass eigenstates. However, this would require the detector to discriminate between neutrinos and antineutrino, which 
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FIG. 8: The neutrino yields as a function of z = E v /m x at the surface of the Earth for annihilation of 250 GeV WIMPs 
into t + t~ in the center of the Earth for the standard oscillation parameters and without neutrino oscillations. The electron 
neutrino yields are not plotted, since they (both with and without oscillations) are identical to the muon neutrino yields without 
neutrino oscillations. 

current neutrino telescope designs cannot do. 

The effects of non-zero #13 and the type of neutrino mass hierarchy have also been studied in Ref. [6| . The neutrino 
flavor ratio results presented in Ref. [6j are in agreement (up to the small effects due to the use of slightly different 
neutrino oscillation parameters) with our low-energy results, where neutrino interactions do not play a major role. 
For the high-energy part of the spectra, it is apparent that such flavor ratios cannot be computed accurately while 
neglecting neutrino interactions. 

VIII. ANNIHILATION IN THE EARTH 

In the case of WIMP annihilations in the Earth, the procedure is the same as for the Sun, i.e., we let the neutrinos 
be produced around the center of the Earth and then propagate them to the detector using the same three-flavor 
framework as used for the Sun although there are a few differences which we have to take into account. First, we must 
consider the size of the annihilation region. This is fairly easy to do given the temperature of the Earth's core and 
we follow the procedure in Ref. [38[ (which is really Eq. ((7|) applied to the Earth). Second, since the density of the 
Earth is much smaller than that of the Sun, B meson interactions are negligible. Third, neutrino interactions on the 
way out to the surface are negligible for the energies of interest here. The geometry is also simpler, since all neutrino 
telescopes are at the same distance from the Earth's core. Hence, the only effect we need to take into account after 
production around the Earth's core are neutrino oscillations on the way to its surface. For each event, we simulate 
the production point and let the neutrino propagate the actual path length to the detector (using the Earth matter 
density profile from Ref. [3611. 

In Fig. [8l we show the yields of neutrinos at the surface of the Earth for 250 GeV WIMPs annihilating into t + t~ 
in the Earth's core. The effects of oscillations are seen as a decrease of the tau neutrino yield and an increase of the 
muon neutrino yield at low energies (below about 50 GeV). Thus, for higher energies, that are of most interest to 
the currently planned neutrino telescopes, neutrino oscillations do not have a large effect on the neutrino yields from 
annihilations of WIMPs inside the Earth. 
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IX. SUMMARY AND CONCLUSIONS 

In this paper, we have studied neutrinos originating from annihilations of WIMPs inside the Sun and the Earth. 
We have computed the yields of all neutrino flavors in a number of different WIMP annihilation channels using an 
event-based full three-flavor Monte Carlo, which we have made publicly available [19]. Our Monte Carlo includes both 
the effects of neutrino interactions and neutrino oscillations. While our study is event-based, earlier studies including 
neutrino oscillations @, S 0] have been performed using a density matrix formalism for the neutrino distributions. 
These two techniques are equivalent, but we have chosen the event-based technique as it is easier to implement as a 
WIMP Monte Carlo for neutrino telescopes. 

We have studied the yields from different annihilation channels for several WIMP masses. However, it would be 
necessary to study the branching ratios and annihilation rates for a specific WIMP model in order to make predictions 
on whether or not that particular model is detectable in a neutrino telescope. In such a study, one would also have to 
consider backgrounds to the WIMP neutrino signal, such as high-energy neutrinos being produced in reactions when 
cosmic rays hit the solar corona (39l . [40L l4lj . The background is of the order of a few neutrino events per year in a 
neutrino telescope such as IceCube. These kind of studies are left for a future work. 

Although our example channel in this paper has generally been the annihilations of WIMPs into pairs of r + and 
t~ for a WIMP mass of 250 GeV, our simulations have been performed for 13 different annihilation channels and 19 
different WIMP masses in the range from 10 GeV to 10 5 GeV. In addition, we have performed all of the simulations for 
different values of the leptonic mixing angle #13, normal and inverted neutrino mass hierarchy, as well as with neutrino 
oscillations turned off. Furthermore, our examples have focused on the results for neutrinos and not antineutrinos, 
which are also accounted for in our Monte Carlo. The full results of our simulations are available at the website in 
Ref. [19(, both as extensive lists of plots of the results, but also as computer- readable data files. In a coming release 
of DarkSUSY [42j , these data files will also be available in a simple-to-use way. 

We have found that, for WIMP annihilations inside the Sun, the main effects of neutrino oscillations arc the 
effective mixing of muon and tau neutrinos during the propagation to the solar surface and the consequent mixing of 
electron neutrinos during the propagation from the solar surface to a distance of 1 AU. However, due to the leptonic 
mixing angle #12 being less than maximal, the mixing between electron neutrinos and the other neutrino flavors is 
not complete. As a result, neutrino oscillations imply that the final yields of muon and tau neutrinos are more or 
less equal at the Earth. This is in contrast to the neutrino yields at production in the Sun, where the electron and 
muon neutrino yields are typically equal, which is then also the case at a distance of 1 AU if neutrino oscillations are 
ignored. 

For non-zero #13, we have found that the adiabatic flavor transitions at the high MSW resonance results in the 
electron neutrinos mixing with the muon and tau neutrinos already during the evolution to the solar surface in the case 
of normal neutrino mass hierarchy. For antineutrinos, we have the opposite situation, since the high MSW resonance 
occurs for antineutrinos only when the neutrino mass hierarchy is inverted. In the case of WIMP annihilations inside 
the Earth, we have found that neutrino oscillations do not play any major role except for at very low energies (less 
than about 50 GeV). 

Our findings imply that, depending on the annihilation branching ratios of the WIMP model under consideration, 
the yield of muon neutrinos, which is the most interesting one from the neutrino telescope point of view, can be either 
increased or decreased. For example, in the case of neutralino dark matter, there is usually a smaller amount of tau 
neutrinos produced. The effective flavor mixing due to oscillations then implies that the muon neutrino yield will 
decrease due to a larger number of muon neutrinos oscillating into tau neutrinos than vice versa. On the other hand, 
dark matter consisting of Kaluza-Klein excitations of Standard Model particles can have large branching ratios into 
charged leptons. Out of these charged leptons, only the taus decay before loosing a significant amount of energy. Thus, 
there will be a larger production of tau neutrinos which can subsequently oscillate into electron or muon neutrinos to 
actually give a larger muon neutrino signal at a neutrino telescope. 
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APPENDIX A: EQUIVALENCE OF MONTE CARLO AND DENSITY MATRIX FORMALISM 



In order for our Monte Carlo simulations to make sense, it is necessary to show that it is equivalent to the density 
matrix formalism for neutrinos from WIMP annihilations presented in Ref. [f|. At any point of propagation r, we can 
construct a density matrix (normalized per annihilation) according to 



(Al) 



where N is the number of annihilations simulated and \v a ) are all of the neutrino states at distance r resulting from 
the Monte Carlo simulations of these annihilations (note that the number of neutrinos does not need to be N). A 
similar definition can be made for the density matrix p(r, E) for antineutrinos. Since the evolution equation is linear 
in the density matrix, it is sufficient to consider the evolution of this density matrix for a pure neutrino state. In this 
case, the density matrix takes the form (up to normalization) 



p(r,E)=S(E-E )p 



(A2) 



where Eq is the energy of the neutrino and po is a matrix in flavor space describing the flavor composition. In our 
Monte Carlo simulations, the same state is described by the state vector \v) such that po = \v) {v\. If disregarding the 
appearance of lower-energy neutrino states due to interactions (although keeping the degradation of the amplitude at 
Eq due to interactions), then the statistical evolution of the state \v) in our Monte Carlo simulations is given by 



dr 



H 



where H is the neutrino flavor evolution Hamiltonian and T — diag(r e 
this follows that the statistical evolution of the density matrix is 



(A3) 

r^jTr) is the interaction rate matrix. From 



^.= (j L H,p ]- 1 -{T,p })ME-E, ) ). 



(A4) 



However, the right-hand side of this equation still lacks the terms from states which have lost energy due to NC 
interactions and from states created due to secondary neutrinos from tau neutrino CC interactions. The effect of both 
of these interactions is to add new states into the collection of states from which we build our density matrix. If the 
NC interaction rate for a neutrino of energy E to produce a neutrino of energy between E and E + dE is j(E , E)dE, 
then the probability of adding a neutrino state in this energy range during the propagation of the small distance Ar 
is Ar"f(Eo, E)dE. In addition, the added state will have the same flavor composition as the original state, i.e., po. 
Thus, it follows that 



dp 
dr 



l(E ,E)p 



(A5) 



NC 



for the case of our single neutrino state. 

For the regeneration of neutrinos with lower energies in CC interactions, we note that neutrinos are only regenerated 
if the interacting neutrino is a tau neutrino. If the interaction rate of tau neutrinos at a neutrino energy of Eq is 
Tq C (Eq), then the rate of tau neutrino interaction for our neutrino is po.rrFJ. c (i?o), since po,rr is the probability for 
our neutrino to currently be in the \v T ) state. If the energy distribution of the new tau neutrino states is f T (Eo, E) 
and the energy distribution of the new \v e>ll ) states is f Btll (Eo,E), then we obtain 



dp 
dr 
dp 
dr 



cc 



cc 



U T po : TTTQ C (E )f T (E , E), 

^-e,ixPO,TT^Cc(Eo)fe,fj.(Eo, E), 



(A6) 
(A7) 



where IL. is a projector onto the \v T ) state and II e M is a projector onto the \v e 
In total, we now have 



subspace. 



dp 
dr 
dp 

dr 



(r,E) = [H,p}~ -{T,p} + ^(E ,E)p a +Tl T p Q , TT r T cc (E )f T (E Q ,E), 
(r,E) = U e ^p ^ T T T cc (E )f e ^(E ,E) 



(A8) 
(A9) 
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for p(r, E) — po5(E — Eq) and p(r, E) = 0. Inserting this into the density matrix evolution equations of Ref. [f|, we 
also obtain the very same evolution equations, which proves that the two treatments arc equivalent on a statistical 
level. 
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